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Abstract 

Accurate prediction of interfacial slip in nanoscale channels is required by many microfluidic applica- 
tions. Existing hydrodynamic solutions based on Maxwellian boundary conditions include an empirical 
parameter that depends on material properties and pore dimensions. This paper presents a derivation of a 
new expression for the slip coefficient that is not based on the assumptions concerning the details of solid- 
fluid collisions and whose parameters are obtainable from equilibrium simulation. The results for the slip 
coefficient and flow rates are in good agreement with non-equilibrium molecular dynamics simulation. 
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Many recent experiments on nanoscale fluid flow of Newtonian liquids 111 |2|] and related sim- 
ulation studies |3, 4,5,^ support Maxwell's prediction of the possibility of molecular slip at a 
gas-solid interface OD. It is known to exist where the Knudsen number, defined as the dimension- 
less ratio of molecular free path to some characteristic length, Kn = A/L, is non-negligible, and 
is usually associated with low densities where A is large. However, in narrow capillaries (where 
L is small), slip can be observed even at liquid densities. In the general case it is characterized 
by a slip coefficient, Zs [tZilSl], which in the absence of temperature gradients relates the collective 
molecular velocity at the wall, the slip velocity, to the shear rate, Ug = LVu DTI]. Using kinetic 
theory Maxwell derived a microscopic expression for the slip coefficient DTD, that can be written as 
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where A = 2r]/pcis the mean free path, r] is the shear viscosity, p is mass density of the fluid, 
and c is the mean speed of the molecules. Considering only two types of wall collision, specular 
and diffuse (Knudsen) reflection, he introduced a coefficient a that defines a fraction of specularly 
reflected molecules. In more broad sense, a defines the fraction of the flux of tangential momen- 
tum transmitted in collisions and is often called the 'accommodation coefficient' (TMAC) [8]. Its 
value is defined by the details of the solid-fluid interactions and ([T]) implies finite slip even for 
purely diffusive reflections (a = 1). For specularly reflecting surface the slip coefficient diverges 
since the fluid cannot grip onto the surface. 



Slip in nanoscale fluid flow depend on many parameters inculding surface roughness [|6D, elec- 
tric properties of the interface [9], wetting conditions js], chemical patterning of the surface [ 101. 



and is a nonlinear function of the dynamic state [lllll . Application of the generalized Navier-Stokes 
(NS) hydrodynamics to problems of fluid flow on the nanoscale is an attractive but highly non- 
trivial task even in simple cases such as plane Poiseuille flow. An accurate solution requires a 
knowledge of the material parameters of the fluid as a function of local density which deviates 
from its bulk value in the proximity of the interface. As a result, the velocity profile in narrow 

r L 

pores also deviates from the macroscopic prediction UM- The usual approach is to regard the 
fluid as incompressible and to replace the complex non-uniform flow problem plus simple no-slip 
boundary conditions with a simple flow problem but with boundary conditions that has been called 
'exceedingly difficult' for theoretical investigation [|81], p. 96. The discontinuity in the flow field 
is introduced in this approach in the same way the surface excess was introduced by Gibbs in 
his treatment of interface boundaries at equilibrium as illustrated in Fig. [T]for the case of plane 



in; 




FIG. 1: Full line, velocity profile for the plane Poiseuille flow with surface slip; dashed line, extrapolated 
velocity profile. Vertical gray line marks the position of the solid. Dash-dotted line denotes the velocity 
gradient at the wall. Also shown are sUp velocity, u^, and slip coefficient, 4- 



Poiseuille flow in z direction (u = {0, 0, u^}) between walls at y = ±h, where the symmetry of 
the solution was taken into account and only one half was plotted. We emphasize that the finite 
slip Ms that appears in this approach is a purely artificial device introduced only to match the ap- 
proximate solution with the solution of full problem in the middle part of the channel. The full 
solution for the continuum velocity field does decay to zero at the wall, as required by continuity 
of stresses in classical hydrodynamics lll3ll . It should be noted however that the limiting value of 
the velocity field cannot be observed in any experiment or in a statistical particle-based simula- 
tion with a continuous solid-fluid potential since the velocity field, which is defined everywhere 
the fluid density is non-zero, cannot be measured from particle velocities in regions where par- 
ticles are not observed due to finite sampling of the statistical ensemble, i. e. where the particle 
Boltzmann factor is diminishingly small but non-zero. 



There were a number of attempts to estimate the TMAC from kinetic theory il4lll5|l and using 



3 



molecular dynamics simulation of simple liquids 



170 . However, the results obtained by 
different groups and using different methods are inconcluisive Jl], with large scatter in obtained 
values of slip coefficient spanning two orders of magnitude. The agreement between recent ex- 



periments with various surfaces [|18 



30,13, 



2111 . and and theory is also poon with experimental 



values for a generally higher than the corresponding theoretical expectations [|15l . 



11. One of 



the limitations of the original Maxwell model is that it was developed to solve the half-space flow, 
or Kramers, problem [8] in the dilute gas regime and faces problems describing flow in pores of 
finite width. It is assumed that the TMAC is a local property independent of the pore width. The 
following argument however demonstrates why it could not be so. Viscous wall stress in plane 
Poiseuille flow that is proportional to the velocity gradient at the wall scales as the pore width 
for the flow of gas of the same density. This stress is transmitted to the wall through collisions 
that are independent of the pore width and therefore the momentum transferred to the wall in each 
collision and consequently the TMAC depends on the pore width. 

In this paper we present a method of calculating the slip coefficient from equilibrium simulation 
that does not require assumption about the type of wall collisions avoiding thus the necessity of 
calculating the TMAC. We illustrate the method on a simple case of gravity driven plane Poiseuille 
flow in a pore of width H = 2h and with acceleration due to an external field g = {0, 0, g} as 
sketched in Fig.[T] The method can be extended to a more general case of Poiseuille flow induced 
by a pressure gradient using the equivalence between the pressure gradient and gravity-driven 
force in the direction of flow [|l3ll . It is convenient to reformulate the slip velocity problem (the 
Dirichlet boundary condition) in terms of interfacial viscosity (Neumann boundary condition), and 
the lateral wall stress. This idea goes back to Navier 1221 ] who obtained the boundary condition for 
the velocity field on the basis of particle arguments (cf. last eqn. on p. 415 in Ref. [22']) as. 
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dy 
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where his parameter (3 is related to the interfacial shear viscosity rf via = r]/ (3. Note that for non- 
linear velocity profiles, as for Poiseuille flow, the slip coefficient, Zj, (shown in figure) is different 
from the slip length, defined as the distance from the wall where the extrapolated velocity profile 
vanishes. However, if the difference between them is small, simple geometric consideration allows 
us to establish the relationship between the two types of boundary condition. 

Neglecting viscous dissipation, the wall shear stress in terms of the external driving force acting 
on fluid particles, which in this case is simply Uy^ = pgh [|23ll . can be equated with the Stokesian 



drag force per unit area exerted on the wall by the moving fluid [12411 . 

Fs Mu _ puh 

where M is the total mass of the fluid in the pore; u = h^^ Uz{y)dy is the mean fluid 
velocity, and the relaxation time r can be calculated in molecular dynamics simulation from 
the Langevin equation for the fluid subsystem considered as a single Brownian particle using 



fluctuation-dissipation theorem. It has been shown recently uM that the fluid velocity autocorre- 
lation function decays exponentially, 

C{t) = M-^kTexp {-t/r) . (4) 

This provides a simple way to determine r in an equilibrium simulation by fitting a one-parameter 
exponential to velocity autocorrelation data. From the two expressions for the wall shear stress, a 
simple relationship between the fluid velocity and the acceleration due to the external force, 

u = rg. (5) 

can be established. This surprisingly simple result shows that within the limits of linear regime 



Ullll the rate of fluid flow in non-equilibrium steady state can be estimated from the characteristic 
time of the decay of fluctuations at equilibrium. Using it, we can also estimate both the slip 
velocity and slip coefficient. Since the hydrodynamic solution in this case is given by a quadratic 
velocity profile with slip, 

uM = ^{h'-y')+u,, (6) 

using Q and the definition of u, one obtains for slip velocity 

n,= (r^^)g, (7) 



37] J 

and by substituting it into the definition of the slip coefficient, = — du/dy\f^Us, one obtains 
finally 

This is the main result. It shows that the slip coefficient is independent of the external force (flux), 
but nonlinearly depends on the pore width, both directly and indirectly through the relaxation time 
T = T{h). The connection with the Maxwell's result can be established by using the relationship 
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between the TMAC and relaxation time, /o, r = (/qo;)"^ MM- Taking the kinetic theory expres- 
sions for the wall collision frequency per particle, /o = c/4/i and viscosity, obtain the expression 
for the slip coefficient 

/s = A- - ^ (9) 

that can be directly compared with Maxwell's result ©. We note here that usage of kinetic theory 
expression for /o introduces a significant error at liquid densities, and it was used in deriving ^ 
only to make the same level of approximation as was in deriving ([T]). 

In order to test the accuracy of dS]), two series of equilibrium molecular dynamics runs were per- 
formed. There are two ways of controlling the Knudsen number in simulation: either by changing 
the mean density of the system (the mean free path) or the pore width (the characteristic length), 
and we used both. In the first series, the relaxation time r was calculated as a function of fluid 
density for a pore of fixed width. In the second, the pore width was varied while the normal wall 
pressure was kept constant. The slip coefficient estimated using dS]) was compared with the values 
obtained directly in the parallel set of nonequilibrium molecular dynamics (NEMD) simulations 
of the Poiseuille flow for the same systems. The system consists of supercritical Lennard- Jones 
(LJ) fluid, confined between two walls modeled by a rigid triangular lattice of atoms situated at 
?/w = ^h, and periodically replicated along x and z axes to avoid edge effects. All interactions 
in the system were of the LJ form, U{r) = As [(a/r)^^ — (cr/r)^] , where e and a are the usual 
energy and length parameters. Their values for the fluid-fluid interactions define the corresponding 
scales, and in the following reduced units JlS] are used, denoted by the astersk. The solid-fluid 
interaction parameters in these units were taken as e^i = 0.4348£: and asf = 0.9462cr, which is 
appropriate for methane gas between carbon 'rare-^s walls'. The surface number density of the 



solid was n^a'^ = 1.105 as in our earlier study t2M where more details about the system and 



the numerical scheme can be found. All molecular dynamics calculations were performed using 



the classical molecular dynamics software package MDL [12611 . The nonequilibrium steady-state 
conditions were realized by placing the fluid in a uniform external field parallel to the walls and 
coupling all fluid degrees of freedom to a Nose-Hoover thermostat at T = 2.026 eA;"^ where k 
is the Boltzmann constant. The simulation cell was of dimensions 20. 18 la x H x L^, where the 
dimension in the flow direction, L^, was scaled with the density to keep the number of particles in 
the system around 2500. The timestep was At* = 7.28 x 10~^ in LJ units and the integration time 
in each case was not less than 3.64 x 10^ (50M steps). The value of the acceleration due to the 
external force was varied between g* = A - 10"^ and g* = 0.04 in order to keep the fluid velocity 



below u* = 0.2. The steady states reported here are known to be well within the Newtonian regime 



[I27I1 and the slip velocity in the linear regime [|ll|] . To simplify the comparison we scaled all flow 
related properties to a common value g* = 4.95 ■ 10^^. 

In the first series three pore widths were considered, H = 5.328cr, 10.499cr, and 20.997cr; (in 
the future denoted as 5a, 10a, and 20a for brevity) and several number densities ranging from 
n* = 0.02 (n* = m^^pa^), corresponding to the rarefied fluid at Kn = 2.34 for a narrow pore 
(where we used the kinetic theory expression for the mean free path and the hard sphere diameter 



for LJ [|28D), to a dense state of n* = 0.8 that for a wide pore gives Kn = 0.015, thus spanning 



more than two orders of magnitude of Knudsen numbers. In all cases the Reynolds number defined 
as the ratio of inertial and viscous forces. Re = puH/r], was kept small. Re < 10. 

In order to calculate slip velocity using ([8]) we need the estimates of the shear viscosity in the 
channel. It can be calculated in equilibrium moleuclar dynamics e.g. from the stress autocorrela- 
tion function using Green-Kubo relations ri29l|. For the bulk LJ fluid it has been recently accurately 
estimated in molecular dynamics simulation [|29l] . To simplify the calculations we used an empiri- 
cal equation of state for the viscosity [|30|] at a density equal to the mean density in the central part 
of the channel, where it is uniform to fit the simulated data of |i29| at the required temperature. 
The accuracy of this procedure was estimated by comparing the obtained values with values cal- 
culated from the NS hydrodynamics using the quadratic fit to NEMD velocity profiles. The results 
for three pore widths are presented in Fig. |2l The agreement for the two wider pores is excelent. 
Deviation from the bulk values for the third pore, H = 5a, is due to the overlap of adsorbed layers 
at two surfaces and as a consequence, to inaccuracy in determination of the corresponding bulk 
density. At low densities the viscosity markedly deviates from bulk values for all three pores when 
Kn > 1. 

The ratio of the slip coefficients to the pore width, calculated in NEMD and estimated from the 
relaxation times and bulk viscosities at densities equal to that in the middle of pore using ^ are 
compared on Fig. |3]for three pore widths. The relaxation time was estimated from the exponential 
fit to the collective velocity autocorrelation function calculated in equilibrium molecular dynamics 
simulaton using the procedure described in uM- The statistical uncertainty in all cases is of the 
order of symbol size and is slightly higher for wider pores since the longer relaxation times in 
this case require more accurate estimation of the velocity autocorrelation time at long correlation 
times. For all pores the slip coeficients calculated using the two routes agree within statistical 
uncertainties. For a given pore, both ([T]) and ^ predict linear dependence of the slip length on 
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FIG. 2: (color online). Comparison of the calculated shear viscosities for three pore widths as a function of 
reduced number density in the centre of pore, symbols, with empirical EOS for bulk fluid, dashed line. 

the Knudsen number with zero and negative offset, correspondingly. For two wider pores the 
behaviour is observed at high Knudsen numbers, Kn > 1, but with positive offset that increases 
with pore width. Surprisingly, for the narrowest pore, the slip length appears to be independent of 
the Knudsen number (density) for Kn > 1. Since density enters only the first term in ([8]), this resit 
indicate that in the low density regime the relxation time in narrow pores increases linearly with 
density. 

The results also show that for the systems with the same Knudsen number the slip length 
increases with the pore width. In order to establish whether there is a limiting value of the slip 
coefficient a second series of calculations was performed at the normal pressure that corresponds 
to reduced density of n* = 0.125 up to the pore width H = 210cr. Obtained values of the reduced 
slip coefficient are compared on Fig. |4] with the Maxwell's results using ([T]). The results estimated 
using ^ agree with those obtained directly in NEMD within statistical uncertainty. They show 
that at about H = 25a the slip coefficient reaches the limiting value 1^ = 4.77(2)(t. Maxwell's 



8 




12 3 

Kn 



FIG. 3: (color online). Slip coefficients calculated directly in NEMD (open symbols and full lines) and 
estimated from the relaxation time (eqn. ([8]l), filled symbols and dashed lines) as a function of Knudsen 
number for three pore widths. Lines are drawn to guide the eye. 

theory, on the contrary, predicts linear scaling of the slip length with pore width and for the widest 
pore studied it gives the value = 77.5a (note the difference in scales). 

This work was funded as part of the National Physical Laboratory's Strategic Research Pro- 
gramme and EPSRC under grant GR/N64809/01. 
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FIG. 4: (color online). Slip coefficients (left scale) calculated directly in NEMD, triangles, and estimated 
from the relaxation time (eqn. ([H)), circles, as a function of pore width. For comparison, slip coefficients 
estimated from the Maxwell's theory, eqn. ([T]), diamonds, are also shown (right scale). Lines are drawn to 
guide the eye. 
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